library(GEOquery)
library(tidyverse)
library(tidySummarizedExperiment)
library(survival)
library(survminer)

slice_expr_quart <- function(df, gene){
  df |>
    mutate(
      goi = {{ gene }},
      expr = case_when(goi <= quantile(goi, c(.25)) ~ 'lowest 25%',
                       goi >= quantile(goi, c(.75)) ~ 'highest 25%',
                       .default = 'NA')) |>
    filter(expr != 'NA')
}

# GSE87371 -------
gset <- getGEO("GSE87371", GSEMatrix =TRUE, getGPL=TRUE, AnnotGPL = TRUE)

gset <- gset[[1]]

gset_meta <- gset@phenoData@data |>
  select(geo_accession, characteristics_ch1.2, characteristics_ch1.3, characteristics_ch1.5, characteristics_ch1.8:characteristics_ch1.12) |>
  mutate(age = str_remove(characteristics_ch1.2, 'age: ') |> as.numeric(),
         sex = str_remove(characteristics_ch1.3, 'Sex: '),
         treatment = str_remove(characteristics_ch1.5, 'treatment: '),
         pfs_time = str_remove(characteristics_ch1.8, 'pfs_time: ') |> as.numeric(),
         cens_pfs = str_remove(characteristics_ch1.9, 'cens_pfs: ') |> as.numeric(),
         os_time = str_remove(characteristics_ch1.10, 'os_time: ') |> as.numeric(),
         cens_os = str_remove(characteristics_ch1.11, 'cens_os: ') |> as.numeric(),
         subtype = str_remove(characteristics_ch1.12, 'coo: '),
         .keep = 'unused') |>
  as_tibble()

probe_gene <- gset@featureData@data |>
  as_tibble()

probe_gene |>
  select(ID,`Gene symbol`,`Nucleotide Title`) |>
  write_csv('DLBCL-Btk/data/GPL570-array-meta.csv')

## also filter by JetSet best probe set
probe_list <- probe_gene |>
  filter(`Gene symbol` %in% c('MYD88','BTK','FCGR2B') | ID %in% c('205297_s_at')) |>
  select(ID, `Gene symbol`)

probe_list <- probe_gene |>
  filter(ID == '225897_at') |>
  select(ID, `Gene symbol`)

### CLCC1 best probe
probe_list <- probe_gene |>
  filter(ID == '213628_at') |>
  select(ID, `Gene symbol`)

tse <- gset |>
  makeSummarizedExperimentFromExpressionSet()

## find high MYD88 expression sample
gene_expr <- tse |>
  assay() |>
  as_tibble(rownames = 'ID') |>
  right_join(probe_list) |>
  pivot_longer(where(is.numeric), names_to = 'sample', values_to = 'expr')

gene_expr |>
  select(-ID) |>
  pivot_wider(names_from = `Gene symbol`, values_from = expr) |>
  left_join(gset_meta, join_by(sample == geo_accession)) |>
  write_csv('DLBCL-Btk/results/GSE87371_DLBCL.csv')

gene_expr <- read_csv('DLBCL-Btk/results/GSE87371_DLBCL.csv')

gene_expr |>
  slice_expr_quart(MARCKS) |>
  survfit(Surv(os_time, !cens_os) ~ expr, data = _) |>
  ggsurvplot(pval = TRUE,
             risk.table = 'nrisk_cumcensor',
             legend.labs = c('Highest 25%','Lowest 25%'),
             title = "GSE87371 DLBCL patients OS with MARCKS expression") +
  xlab('Survival time (month)')

gene_expr |>
  slice_expr_quart(MARCKS) |>
  survfit(Surv(pfs_time, !cens_pfs) ~ expr, data = _) |>
  ggsurvplot(pval = TRUE,
             risk.table = 'nrisk_cumcensor',
             legend.labs = c('Highest 25%','Lowest 25%'),
             title = "GSE87371 DLBCL patients PFS with MARCKS expression") +
  xlab('Survival time (month)')

## only R-CHOP treated patients
gene_expr |> filter(str_detect(treatment, 'CHOP')) |>
  slice_expr_quart(MARCKS) |>
  survfit(Surv(os_time, !cens_os) ~ expr, data = _) |>
  ggsurvplot(pval = TRUE,
             risk.table = 'nrisk_cumcensor',
             legend.labs = c('High MARCKS','Low MARCKS'),
             title = "GSE87371 DLBCL patients OS treated with R-CHOP") +
  xlab('Survival time (month)')

gene_expr |> filter(!str_detect(treatment, 'CHOP')) |>
  slice_expr_quart(MARCKS) |>
  survfit(Surv(os_time, !cens_os) ~ expr, data = _) |>
  ggsurvplot(pval = TRUE,
             risk.table = 'nrisk_cumcensor',
             legend.labs = c('High MARCKS','Low MARCKS'),
             title = "GSE87371 DLBCL patients OS with other treatment") +
  xlab('Survival time (month)')

## save tidy data ----------
probe_gene <- probe_gene |>
  select(ID, `Gene symbol`)

tse |>
  assay() |>
  as_tibble(rownames = 'ID') |>
  left_join(probe_gene) |>
  relocate(ID, `Gene symbol`) |>
  write_csv('mission/DLBCL-Btk/data/tidy/GSE87371_expression.csv')

gset_meta |>
  write_csv('mission/DLBCL-Btk/data/tidy/GSE87371_metadata.csv')

# GSE32918 -----------
gset <- getGEO("GSE32918", GSEMatrix =TRUE, getGPL = TRUE, AnnotGPL = TRUE)

gset <- gset[[1]]

gset_meta <- gset@phenoData@data |>
  as_tibble() |>
  select(title, geo_accession, characteristics_ch1.1, characteristics_ch1.3, characteristics_ch1.4, characteristics_ch1.6, characteristics_ch1.7) |>
  mutate(age = str_remove(characteristics_ch1.3, 'age: ') |> as.numeric(),
         vital_status = ifelse(str_detect(characteristics_ch1.6, 'Alive'), 0, 1),
         subtype = str_remove(characteristics_ch1.1, 'predicted class: '),
         sex = str_remove(characteristics_ch1.4, 'Sex: '),
         years = str_remove(characteristics_ch1.7, 'follow-up years: ') |> as.numeric(),
         title = str_remove(title, '_.+'),
         .keep = 'unused')

ambi_sample <- gset_meta |>
  select(title, geo_accession)

gset_meta <- gset_meta |>
  distinct(title, .keep_all = TRUE)

## array annotation ------
probe_gene <- gset@featureData@data |>
  as_tibble()

probe_list <- probe_gene |>
  select(ID, Symbol) |>
  filter(Symbol %in% c('CLCC1')) |>
  distinct(Symbol, .keep_all = TRUE)

## find expression matrix of interest
tse <- gset |>
  makeSummarizedExperimentFromExpressionSet()

myd_expr <- tse |>
  assay() |>
  as_tibble(rownames = 'ID') |>
  right_join(probe_list) |>
  pivot_longer(where(is.numeric), names_to = 'sample', values_to = 'expr')

myd_expr |>
  left_join(ambi_sample, join_by(sample == geo_accession)) |>
  summarise(expr = mean(expr), .by = c(Symbol, title)) |>
  pivot_wider(names_from = Symbol, values_from = expr) |>
  left_join(gset_meta) |>
  write_csv('DLBCL-Btk/results/GSE32918_DLBCL.csv')

myd_expr <- read_csv('DLBCL-Btk/results/GSE32918_DLBCL.csv')

myd_expr |>
  slice_expr_quart(MARCKS) |>
  survfit(Surv(years, vital_status) ~ expr, data = _) |>
  ggsurvplot(pval = TRUE,
             risk.table = 'nrisk_cumcensor',
             legend.labs = c('Highest 25%','Lowest 25%'),
             title = "GSE32918 DLBCL patients OS with MARCKS expression") +
  xlab('Survival time (years)')

low_fcgr |>
  survfit(Surv(years, vital_status) ~ FCGR2B_expr, data = _) |>
  ggsurvplot(pval = TRUE,
             risk.table = 'nrisk_cumcensor',
             title = "GSE32918 DLBCL patients OS with FCGR2B expression")

gset_meta |>
  write_csv('mission/DLBCL-Btk/data/tidy/GSE32918_metadata.csv')

probe_gene <- probe_gene |>
  select(ID, ILMN_Gene)

tse |>
  assay() |>
  as_tibble(rownames = 'ID') |>
  left_join(probe_gene) |>
  relocate(ID, ILMN_Gene) |>
  write_csv('mission/DLBCL-Btk/data/tidy/GSE32918_expression.csv')

# GSE31312 ----------
gset <- getGEO("GSE31312",
               GSEMatrix = TRUE, getGPL = TRUE, AnnotGPL = TRUE)

gset <- gset[[1]]

write_rds(gset, 'mission/DLBCL-Btk/data/gse31312.rds')

gset <- read_rds('mission/DLBCL-Btk/data/gse31312.rds')

## array annotation ------
probe_gene <- gset@featureData@data |>
  as_tibble() |>
  select(ID, `Gene symbol`)

## find expression matrix of interest
tse <- gset |>
  makeSummarizedExperimentFromExpressionSet()

tse |>
  assay() |>
  as_tibble(rownames = 'ID') |>
  right_join(probe_gene) |>
  relocate(`Gene symbol`) |>
  write_csv('mission/DLBCL-Btk/data/tidy/GSE31312_expression.csv')
